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Abstract 

We develop a framework for spectrum sensing in cooperative amplify-and-forward cognitive radio 
networks. We consider a stochastic model where relays are assigned in cognitive radio networks to 
transmit the primary user's signal to a cognitive Secondary Base Station (SBS). We develop the Bayesian 
' optimal decision rule under various scenarios of Channel State Information (CSI) varying from perfect 

' to imperfect CSI. In order to obtain the optimal decision rule based on a Likelihood Ratio Test (LRT), 

£N| , the marginal likelihood under each hypothesis relating to presence or absence of transmission needs 

to be evaluated pointwise. However, in some cases the evaluation of the LRT can not be performed 
analytically due to the intractability of the multi-dimensional integrals involved. In other cases, the 
distribution of the test statistic can not be obtained exactly. To circumvent these difficulties we design 
two algorithms to approximate the marginal likelihood, and obtain the decision rule. The first is based 
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on Gaussian Approximation where we quantify the accuracy of the approximation via a multivariate 
version of the Berry-Esseen bound. The second algorithm is based on Laplace approximation for the 
marginal likelihood, which results in a non-convex optimisation problem which is solved efficiently via 
Bayesian Expectation-Maximisation method. We also utilise a Laguerre series expansion to approximate 
the distribution of the test statistic in cases where its distribution can not be derived exactly. Performance 
is evaluated via analytic bounds and compared to numerical simulations. 

Index Terms 

Cognitive radio, Cooperative spectrum sensing, Likelihood Ratio Test, Laplace method, Laguerre 
polynomial, Berry-Esseen theorem, Bayesian Expectation Maximization. 
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I. Introduction 

In recent years, cognitive radio [1], [2] has attracted intensive research due to the pressing demand 
for efficient frequency spectrum usage. In a cognitive radio system, the secondary users (SU) try to find 
"blank spaces", in which the licensed frequency band is not being used by the Primary Base Station 
(PBS). A key requirement in cognitive radio is that the SUs need to vacate the frequency band as quickly 
as possible if the corresponding Primary User (PU) emerges. 

Spectrum sensing is a mandatory functionality in any CR-based wireless system that shares spectrum 
bands with primary services, such as the IEEE 802.22 standard [3]. This standard proposes to reuse 
vacant spectrum in the TV broadcast bands. There has been a significant amount of research on spectrum 
sensing for cognitive radio, see [4], [5] for overviews. 

Essentially, spectrum sensing can be cast as a decision making or classification problem. The secondary 
network needs to make a decision between the two possible hypotheses given an observation vector: that 
the frequency band is either occupied or vacant. The more knowledge we have on the nature of the primary 
user's signal, the more reliable our decision. If no knowledge is assumed regarding the primary user, 
energy detector based approaches (also called radiometry) are the most common way of spectrum sensing 
because of their low computational complexity. Cooperative networks can improve the performance of the 
network by enabling users to share information and create diversity. This helps to combat the detrimental 
effect caused by the fading channels. In this context, cooperative spectrum sensing has been studied 
extensively as a promising alternative to improve the sensing performance. In [6], the authors proposed 
algorithms to optimise detection performance by operating over a linear combination of local test statistics 
from individual secondary users. In [7], the performance of cooperative spectrum sensing was derived. It 
was found that the optimal decision fusion rule to minimize the total error probability is the half-voting 
rule. In [8], centralized and decentralized detection schemes were developed. 

In contrast to those methods, our system model for cooperative spectrum sensing contains the practical 
scenario of channel uncertainty. This includes the case of partial CSI knowledge at the SBS or the more 
severe case of blind spectrum sensing. We also assume that the relays have no processing capability, 
therefore are not capable of performing any local decisions. This is a practical scenario encountered in 
many relay networks [9], [10]. In order to perform LRT, the SBS performs a hypothesis test to decide 
whether the PBS is transmitting or idle in a given frame. As we show, the densities involved in making 
a decision under this framework are intractable, meaning they can not be evaluated point wise. This is 
due to the fact that they involve multi-variate integrals which can not be solved analytically. 
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Contribution: 

1) We propose a novel statistical model to address the problem of spectrum sensing with partial CSI. 
To the best of our knowledge, a cooperative spectrum sensing, where both the channels from the 
PBS to the relays and the channel between the relays to the SBS are only partially known has 
not been addressed previously. In most cooperative CR systems, the relays perform a local soft or 
hard decision and then report their summary statistics to the SBS [11], [12]. In the system model 
we present all the statistical processing is performed at the SBS, thus removing the computational 
complexity from the relays and placing it at the SBS, enabling the use of standard relays systems 
already developed and in operation, making such an approach widely applicable . 

2) We derive the probabilities of detection and mis-detection as well as the associated optimal tests 
under several different scenarios. Some bounds have exact closed-form expressions while others 
have closed-form approximations that we derive via Laguerre series expansion. 

3) For the most complicated case of imperfect CSI, we derive two low complexity algorithms to 
perform spectrum sensing: 

i. The first is based on Gaussian approximation via moment matching. This results in a simple 
closed-form test statistic, for the decision process. In addition we study the approximation 
error providing closed form expression for the bound via Berry-Esseen theorem. 

ii. The second is based on the Laplace approximation of the marginal likelihood which involves 
solving a non-convex optimisation problem. 

The paper is structured as follows: in Section II the stochastic system model is developed and the 
Bayesian estimation problem is presented. Section III presents an analysis of the case of perfect CSI. 
In Section IV we develop the optimal decision rule and approximate the performance for the case of 
imperfect-perfect CSI case. In Section V we present two novel algorithms to perform the hypothesis test 
in the case of imperfect-imperfect CSI. Section VI presents extensive simulation results. Conclusions are 
provided in Section VII. 

Notation. The following notation is used throughout: random variables are denoted by upper case 
letters and their realizations by lower case letters, and bold case will be used to denote a vector or matrix 
quantity. 

II. Problem definition and System model 

In general, in cooperative spectrum sensing model, the task of the SBS is to discriminate between two 
hypotheses, the null (Ho) that the bandwidth is idle versus the alternative (Hi) that the bandwidth is 
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occupied, given a set of observations. Based on the decision regarding the presence or absence of primary 
user's activities, the cognitive radio can utilise the spectrum or vacate it. Since we will formulate the 
problem of deciding whether the channel is occupied or not based on the observational evidence via a 
nested model structure, we can consider the likelihood ratio test framework, see [13]. 

In this paper, the challenging aspect of this problem, that extends beyond solutions developed pre- 
viously, is that the probability distribution (PDF) of the test statistic that we wish to use for inference 
is not known in closed form. In general, it will also depend on a set of unknown parameters, for each 
hypothesis. Therefore we resort to formulating analytic approximate solutions for the distribution of the 

LRT under both hypotheses in order to perform inference. 
A. Statistical model 

The network architecture we consider is a centralized network entity such as a base-station in infrastructure- 
based networks (see Fig. 1). We consider a frame by frame scenario where one PBS may be active 
(transmitting data) or idle (not transmitting) during a frame. If active, its signal is transmitted over 
independent wireless channels and is captured by M relay links. Each relay, instead of making individual 
decisions about the presence of the primary user, simply transmits the noisy received signal to the SBS 
over a fading channel. The SBS is equipped with N receive antennas. We further assume that the SBS 
has only limited knowledge of the CSI (noisy channel estimates), which is a practical scenario [9]. We 
now outline the system model and associated assumptions. 
Model Assumptions: 

1) Assume a wireiess network with one PBS equipped with a single antenna, that may be active or idle 
in a given frame. 

2) In case that the PBS is active, it periodically transmits pilot signals, s(l), I = 1, . . . ,L, within a 
frame of L symbols, see [3]. This model assumption will be discussed in Remark 2 below. 

3) At each frame the received signal at the m-th relay ( m = 1, . . . , M) is a random variable given as a 
composite model, where T-Lq and %\ correspond to idle and active model hypotheses, respectively: 



where F m (l) denotes the channel coefficient between the PBS and the m-th relay and V m {l) is the 
unknown noise realization associated with the m-th relay receiver. Note, each of the relays is equipped 
with single receive and single transmit antenna. 
4) The relays re-transmit their received signal, {R m (l)} m=1 , over M fading channels. These channels 
can either occupy the same frequencies as the PBS-Relay channels or be dedicated reporting channels. 
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5) The received signal at the SBS from all M relays at epoch Z can be 



written as 
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where Y(Z) E c^xi 2 S j/j e received signal at the l-th sample, G(Z) G C NxM is the random channel 
matrix between therelay and the SBS, F(Z) = [Fi(l), ■ ■ ■ ,F M (l)] T £ C Mxl is the random channel 
vector between the PBS and the relays. The random vector, W(Z) £ C Nxl , is the random additive 
noise at the SBS, and V '(I) = [Vi(Z),--- ,V M {l)] T G C Mxl is the random additive noise at the relays. 
Remark 1: Note, the dimension N at the SBS can be attributed to several factors. For example, the SBS 
can be a MIMO receiver equipped with N receive antennas. A different option is that the relays, while 
observing the same frequency band, transmit their information over M dedicated orthogonal frequency 
bands (i.e. N = M). In that case, G would be a diagonal matrix. Here, we wish to make the system 
model as general as possible and not impose particular constraints or assumptions. 
Remark 2: Cognitive radio standard defined in IEEE 802.22 is implemented in the TV bands [3 J. The 
TV bands digital signals can be either ATSC (North America), DVB-T (Europe), or ISDB (Japan). These 
standards contain within them many features, such as pilot symbols and synchronization patterns. For 
example, ATSC signals [14] contain a 511-symbol long PN sequence, pilot symbols and synchronization 
patterns of 828 symbols. This makes our assumption regarding pilot symbols and synchronized transmis- 
sion practical. 

Remark 3: CSI can be obtained using the knowledge of the pilot symbols from Remark 2. If the relays 
have the capability of performing channel estimation, they can forward these estimates to the SBS. The 
SBS can also perform channel estimation to obtain matrix G. 
Prior specification: 

Here we present the relevant aspects of the Bayesian model and associated assumptions. 
1) The PBS is active or idle with prior probabilities P (Hi) and P (Ho), respectively. 
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2) All the channels are time varying, meaning that they are constant within a symbol, but may change 
from one symbol to the next. 

3) The SBS has only a noisy estimate of the true channel realisation, G. This is the result of a channel 
estimation phase which we do not detail here, see [15]. A common approach is to model the channel as 
G(Z) = G(Z) + A, where G(Z) is the noisy channel estimate, and A is the associated estimation error. 
The distribution of G(Z) conditioned on G(Z) and A can be written as G(Z) ~ CN ( G(Z), Eg) , £g 
is the covariance matrix with known elements c G , see details on noisy channel models in [15], [16]. 

4) The SBS has only a noisy channel estimate of the true channel realisation, F(Z). As with the G 
channels, F(Z) can be written as F(Z) ~ CN ( F(Z), Sf) , where F(Z) is the estimated channel and 
Ep = erpl is the covariance matrix with known elements up. Note: this stochastic model covers the 
case where CSI is unavailable, and only the channel prior distributions are available. In that case, 
F(Z) = and of = 1. 

5) The additive noise at the relays is a zero-mean i.i.d Complex Gaussian distribution, V(Z) ~ CN (0, Sy) , 
where Sv = o v I is the covariance matrix, known at the SBS. 

6) The additive noise at the SBS is a zero-mean i.i.d Complex Gaussian distribution, W(Z) ~ CN ( 0, £w) , 
where Ew = °wl 1S covariance matrix, known at the SBS. 

7) The symbols s(Z) are known at the SBS in the form of pilot symbols [3]. For ease of presentation and 
w.l.o.g, we assume that s(l) = 1, VZ G {1, . . . , L}. 



B. Spectrum sensing decision criterion 

The objective of spectrum sensing is to make a decision whether the spectrum band is idle or active 
(choose Hq or 1-Li) in a given frame, based on the received signal at the SBS. To solve the decision 
problem, we will take a Bayesian approach. That is we consider the B ayes' risk formulation of the 
decision problem which generalises the LRT to a Bayesian framework. The problem of designing the 
decision rule can be treated as an optimization problem whose objective is to minimize the cost function: 



(4) 



C = P (Ho) ( Coo [ p(yi:L\n )dy 1:L + C 10 [ P (y 1 ;L\Uo)dy 1 
+ P(Hx)(c if p(yi:L|«i)^i:L + Cn f p (yi :i |«i) dy 1:L Y 

\ JAo JAi J 

It can be shown, [13], that the optimum decision rule is a likelihood-ratio test (LRT) given by 

hi 

Arv , A p(yi#i) > P(Ho)C 10 -Coo A 

{ :L) P(yi-.L\H ) < P(Hi)C 01 -C u 7 ' W 

/L0 
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where C xy is the predefined associated cost of making a decision T~L X , given that the true hypothesis is 
M y , and we define random matrix Yi-x = [Y (1) , . . . , Y (L)]. 

Under both hypotheses Hq and H\, all random quantities in the model are independent. We can 
therefore decompose the full marginals under each hypothesis, piyi-.i^Hk), fe = 0, 1, as 

L 

p(yi:L\n k ) = Y[p(y(i)\n k ). (6) 

1=1 

This decomposition is useful as it allows us to work on a lower dimensional space, resulting in efficiency 
gains for the algorithms we develop in the next sections and requiring no memory storage for data. 

Table I presents a summary of the different scenarios that will be covered in the next Sections as well 
as the type of solution that is provided under each scenario. 

III. Perfect Knowledge of PBS-relays and relays-SBS Channels 

We consider the situation of perfect CSI of both G(l) and F(Z) for all I, which corresponds to Case 
I in Table I. We derive the optimal decision rule and the probabilities of detection and false alarm. 
This scenario enables us to obtain a lower bound on the overall system performance in terms of error 
probabilities in analytic form. 

Lemma 1: The marginal likelihood under perfect CSI is: 



|g(0,f(/)~F(y(0|g(/),f(0) 



civ(o, £(/)), n 

(V) 

CN(fi(l),X(l)),Hi, 



where E(l) 4 <%g(l)g(l) H + a^I, and M (0 = g(Qf (0- 

Next, using the decomposition property of (6), the test statistic and associated decision rule are presented. 
Theorem 1: Under perfect CSI, the optimal decision rule defined in (5) is given by 

a /y > 4 Piyiam ex P -i^ 1 (yw-MO)^- 1 w(y(o-M(0) 

1 1:L) p{yi:L\n Q ) exp-iSf =1 y(0»s-(0y(0 ' ( j 

which results in the following decision rule 



L n ° 



^ 1=1 Til 1=1 



(9) 



where we identify the test statistics according to 

L 

T(Y 1:L ) 4 Y,Re [n(l) H Tr\l)Y{l)] , (10) 



i=i 
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and the threshold for the critical region is given by 

1 L 



(11) 



1=1 



Proof: Using Lemma 1 and the definition of the LRT it follows in the log domain that 

L L 

log A (Y 1:L ) = - Y(Z)^E- 1 (Z)Y(/) - - J2 (Y(0 - n(l)f ^(O (Y(l) - M (J)) 



l=i 



i=i 



L 1 
^Re [ / ,(/)^S- 1 (/)Y(/)] - -^E-^ZKZ). 



(12) 



l=l 



This result is useful as it will provide a lower bound for the achievable Type I (false detection) and Type 
II (false alarm) probabilities as a function of SNR. It can therefore be used as a comparison for our 
approximation results when only partial CSI is known. 

Theorem 2: The probability of detection and false alarm under perfect CSI are expressed analytically 



as 



Pd = p(T(Y 1:L )>T\ni) = Q 



'Ef=iM0"E-i(0M0 



(13) 



and 



P f =p(T(Y 1:L )>T\n ) = Q 



(14) 



respectively, where Q \x] = -A= f°° exp - ' 2 / 2 dt. 

V 2n J x 

Proof: To obtain this we derive the distribution of the test statistic utilised in the Bayes risk criterion 
under both the null and alternative hypotheses as follows 

T(Y 1:L )|W ~ N ( 0, i^Re [(^Stf)" 1 ) £(/) (MO*^) -1 )*] ) = TV ( 0, i £ MO^CO'VO ) , 



i=i 



i=i 



T(Y 1:L )\H 1 ~ JV ( £>e [ M (0 H E(0"V(0] >^X> [(MO^I)" 1 ) E(Q (MO^G) -1 )* 
Vj=i i=i 

L 



n £ m(0 h e(0" V(0, J E M0"e(0~ V(0 



i=l 



Next we establish that adding receive antennas translates into a better overall detection performance. 

Theorem 3: Under perfect CSI, the probability of false alarm, pt, and misdetection, 1 — pj, can be 
shown to be monotonically strictly decreasing as the number of SBS antennas N increases. 
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Proof: To obtain this result, we apply Theorem 1 and show that pd(N + 1) < pd(N). Consequently, 
we prove that fi^E]^ pm < ^n+i^n+i^n+i, where the subscript N, (N + 1) refers to the number of 
receive antennas. We begin by proving that 

SN {vvgNgN + O-W 1 ^) 1 SJV -< g^+1 {^VgN+lgN+l + OwIiV+l) 1 SN+1- (15) 

Consider the following linear model 

Z N = EN X + W, (16) 

where W ~ CN (0, g^In), X ~ CN (0, o\1k), and known mixing matrix g N G C NxK . Then, using 
the properties of the Bayesian MMSE, its MSE covariance matrix, cn, can be written as 

c N = a^I K - dvg^ (ovgivg/v + crw 1 ^) 1 EN- (17) 

If we had an augmented model with (N + l) observations, i.e. Z N G C (Ar+1)xl , g^ G C (JV+1)xM , W N G 
c (N+i)xi t then 

cat+i = g^Ir - (ovgw+igw+i + o-w^+i) 1 gJV+i- (18) 

It is well known that the MSE covariance is strictly decreasing with the number of observations [13], 
that is 

x H gf (ovgjvgw + ow 1 ^) 1 gTvx < x H gN +1 (civgTV+ig^+i + a^I N+1 ) 1 gw+ix, (19) 

and in particular, consider x = f G C Mxl . ■ 
Remark 5: Theorem 3 does not hold for the number of relays, i.e. it's not necessarily true that for 
fixed N and L, p d (M + 1) < p d (M). 

IV. Imperfect PBS-relays CSI and Perfect relays-SBS CSI 

In this section we consider the system model under which we can assume that the receiver has perfect 
CSI of G(l) but only partial knowledge of F(Z), which corresponds to Case II in Table I. We derive the 
optimal decision rule and the probabilities of detection and false alarm via Laguerre series expansion. 

Lemma 2: The marginal likelihood under Imperfect PBS-relays CSI and Perfect relays-SBS CSI is: 

CN (o,s % m), H 
Y(0|g(0~*Xy(0lg(0)= (20) 

[ CN(n(l),X ni (l)),Hi, 

where E^(Z) A o^g(l)g(l) H + a^I , E Hl (l) 4 + 4) g(Z)g(0 H + a^I 
and n(l) = g(/)F(Z). 
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Having stated the likelihood in Lemma 2 we then present the corresponding test statistic for the Imperfect 
PBS-relays CSI combined with perfect relays-SBS CSI setting. 

Theorem 4: Under Imperfect PBS-relays CSI and Perfect relays-SBS CSI, the optimal decision rule 
in (5) is given by 

i 



A(Yi :L ) 



tt l i , rz ~-iMi)-rfi)) H vz\(i)(y(i)-rfi)) 

p(yi:L|fti) nf=iP(y(Q|fti) Ui=1 M w lE(i )wi |^ exp ' 

p(yi:L\n ) nf=iP(y(OI«o) " 



1A '=i (2.)^|e„ (o| 1/2 exp 



which results in the following decision rule 
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+ ^(a(0"a(0+M0^%(0M0) < ^|c(0Y(0 + a(0| 2 . 

T^i z=i 



//ere we identify the test statistics according to 



T(Y 1:L ) = J>0)Y(0 + a(0| 2 . 



(21) 



z=i 

a«<i f/ie resulting Bayes risk threshold is defined as 



r = log 7 + ^2 lo s 



Z=l 



+ X)(a(I) flr a(Z) + /i(0^;(0M(0), 



vwY/i c(l) H c(l) = (E^(Z) - S«}(0). a W - c(l) -1 £«}(0M0- 

Proof: Using the result in Lemma 2 and the definition of the LRT, produces 



21ogA(Y 1:i ) = ]T>g 
i=i 

L 

= ]T log 



*W0 



z=i 



WO 
WO 



WO 



£ Y(i)*£^ (l)Y(O - <T (Y(0 - (I) (Y(0 - M (0) 

(=1 1=1 

L 

|c(0Y(i) + a(/)| 2 - a(0 H a(0 - M0^ (z>(0- 



Z=l 



The distribution of the test statistic in (21) is asymptotically % 2 in L. However, in practical systems the 
number of frames is typically small and therefore this asymptotic result can not be applied. Therefore, 
in these cases the distribution of the test statistic in (21) is not attainable in closed form and deriving 
the probability of detection and false alarm needs to be approximated. The reason for this is that the 
L-fold convolution of non-central \ 2 random variables, each with different centrality parameter, can not 
be solved in closed form. 

There does however exist a rich statistical literature approximating the distribution of the linear 
combination of non-central \ 2 random variables. The solutions to finding an approximation to the PDF 
and CDF of a linear combination of non-central \ 2 involve a range of series expansions, saddle point 
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approximation type methods and the Weiner Germ modifications, see in depth discussions in [17] and 
[18]. In this paper we consider the Laguerre series expansion distributional approximations to attain 
the probability of false alarm and mis-detection. This class of approximation has tight error bounds 
represented as a function of the order of the series, see [19]. 

The Laguerre series expansion for the probability of false alarm and mis-detection is characterized by 
the parameters: p the order of the series expansion; /j,q a parameter that controls the rate of convergence 
of the series expansion; and /3 values selected to control the error of approximation for a given p, see 
discussion in [19]. In addition, this approximation has the property that for different settings of the 
parameter /xq we can obtain other series expansions in the literature such as setting n$ = v/2 = p which 
gives the expansion of [18]. 

Theorem 5: Under Imperfect PBS-relays CSI and Perfect relay s-SBS CSI, the probability of false 
alarm and miss-detection are approximated by the generalized Laguerre series as: 

p d = p (T(Y 1:L ) > T\Hi) = 1 - F p (r|Hx) (22a) 
p f = p (T(Y 1:L ) > r|7^o) = 1 - F P (T\H ), (22b) 
where V is the Bayes risk threshold given in Theorem 4 and 

x™ ~ . x™> - ( 4l, r( ;;; +1) g^4^ (^) . ^ > »,< EK , 

(23) 

with coefficients having the recurrence relations in the setting //q > and p = v/2 + 1 given by 



,y W2+1 f_lv-£ sm^p-jtp) \ 

m = 2(- + l exp( !Ltl »™+^"o)j!: TT (/3 Mo + a Ap - a )) u(l)l . 



j fc-i 

m k = - ^2 mjdk-j k > 1, 

3=0 

2/Uo \PMo + ai(p- po)J 2 \ (j[i Q + a t (p - p J 



The corresponding PDF is given by 
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with p = u/2, and v = Y2i=i ^(0 an d the following recurrence relations for the coefficients, 

f/2 / , m , \ L / \ -MO/2 



co 



fe-i 



i a k—ji 



3=0 



+Ef (i 



1 - a ; /3 



J > 1, 



+ (anP)(p/iio - 1) 

ant/ Ly'(t) = Y^m=o Cj-m - m\ > is the generalized Laguerre polynomial. In addition the 

generalized Laguerre polynomials can be obtained by recurrence relationships, 

jLf\t) = (2j + a - 1 - t) 41 (t) - (j + a - 1) L<f] 2 (t), 

lS(*) = 0, 4 o) (*) = l. 
Proof: To derive this result involves consideration of the distribution for a linear combination of 
non-central x 2 random variables. Consider the identity for LRT statistic given in Theorem 5 as 



1=1 



with 



T(Y 1:i ) = J]||c(0Y(0+a(/)|| 2 = ^||Y 
i=i i 

( \ 



(25) 



Y(l)\n ~ CN 



Y(0|^i ~ CN 



a(0,c(0E Wo c(0 



c(/H0 + a(/),c(0S Wl c(0 



V 



(26a) 



(26b) 



-Hi 



/ 



To obtain the distributional approximations, we require a linear combination of independent x 2 non-central 
random variables. To achieve this for each symbol we apply the following rotational transformation, 
based on SVD decomposition of £w*(Q = U(l)K(l)U T (I) giving transformed random vectors with i.i.d 
elements Z(l) = U(Z)A~ 1/2 (0Y(Z) ~ CN(U(l)A- l / 2 (l)a(l),I). As a result, we obtain a univariate 
linear combination of squared Gaussian random variables, 



LN 



T(Z 1:L ) = 5>iZ 2 (Z), 



(27) 



with ai > a positive weight. Each resulting independent scalar random variable Z 2 (l) ~ X^,(5(l)) w i tri 
non-centrality parameters 5(1). 



January 19, 2013 



DRAFT 



13 



Therefore, under the transformed observation vectors Zi : l, one can obtain the distributions of the test 
statistic in Theorem 5. The equivalent Bayes risk threshold for the transformed data can be easily obtained 
by replacing c(Z) with c(l) = U(/)A~ 1 / 2 (/)c(/) and replacing a(Z) with a(Z) = U(/)A~ 1 / 2 (/)a(/). ■ 



The result of Theorem 5 provides the means to approximate the critical region of the decision rule 
for any observed test statistics for any number of frames. This means we can quantify the mis-detection 
and false alarms rates analytically as a function of the number of frames and the SNR with known error 
bounds on the order of approximation. 

V. Imperfect PBS-relays CSI and Imperfect relays-SBS CSI 

In this section we consider the case where the SBS has only partial CSI of both G(7) and F(Z), which 
corresponds to Case III and Case IV in Table I. 

We first consider two scenarios which have practical interpretations before moving to the more general 
case. The first involves consideration of line-of-sight transmissions; and the second assumes high SNR 
scenario, which for both we obtain analytic expression for the marginal likelihood and therefore an 
analytic for the LRT in (5). 

For the ra-th element in (2), after omitting the time dependency I, we obtain: 



M 



y-(n) _ Q(n,m)-p>{m) _|_ y^{n) ^ 



(28) 



m=l 



The m-th term in the summation above can be expressed as: 



Q{n,m) R {r, 



K 



K 



^f(n,m) 



R 



n 



^f(n,m) 



R {m) 



G (n, 



I G {n ' m) 1 R (m) 
(29) 

Each of the terms in (29) form a product of independent Normal random variables. The distribution of 
this product was first derived by [20] and later studied by [21] and the resulting density and Moment 
Generating Function (MGF) are given as follows. 

Lemma 3: The distribution of a product of two independent normally distributed variates Z = XY, 
where X ~ N (X, a^) and Y ~ N (Y, ay) is the solution of the following integral 

1 



p(z) 



oo POO 



exp 2 "x exp 2 °y S (z — xy) dxdy. 



OO J — oo 



llTO'xO'Y 

The Moment Generating Function (MGF) of Z can be expressed as [21] 



M z (t) 



exp 



{p 2 x +pj-)t 2 +2p x p Y t ~ 
2(l~t 2 ) 



(30) 



(31) 
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where px = and py = 

For the special case where X = Y = the integral can be solved analytically as [21 ] 



p(z) 



(32) 



irO~x°~Y 

where Kq (■) is the modified Bessel function of the second kind. 

Following this result of the MGF we obtain an asymptotic result for the marginal likelihood in (6) 
under both hypotheses. 

Theorem 6: The distribution ofY^ in (28) is asymptotically Normal when either pq = > oo 

or pp = > oo. Therefore 



CN ( mo ,Z Ho ) ,H 

(33) 

CN ( mi ,Z Hl ) ,Hi. 



y(n) = G( n ' m )i?( m ) + ~ < 

m=l 

Therefore, assuming only the linear dependence between the N components of Y and ignoring any tail 
dependence in the joint multivariate distribution, we conclude that Y is multivariate Gaussian. 

Proof: See Appendix. ■ 
Theorem 6 shows that under the following conditions, the Gaussian Approximation (GA), that will be 
presented next, is valid: 

1) The CSI estimation error, quantified by ctq and/or dp is low. 

/ j~ YYl ) ( 111/ j 

2) The mean value of one or both of the channels estimates (G ' , R ) is large, i.e. strong line- 
of-sight, for example in Rician channels. 

Next we consider generalising the analysis to relax the assumptions in Theorem 6 making the resulting 
solution widely applicable. Consequently, the distribution of the marginal likelihood in (6) under both 
hypotheses is intractable. This is because it involves finding the distribution of in (28) which can 
not be obtained analytically. This is due to the fact that (Mz (t)) M ^ M^ M z m (*)> which means that 
this distribution is not closed under convolution. 



A. Gaussian Approximation via Moment Matching 

We derive a low-complexity detection algorithm that is based on moment matching so that the dis- 
tribution of the received signal is approximated by a matrix variate Gaussian distribution based on the 
results obtained in Theorem 6. We show under which conditions this approximation is valid and asses 
the approximation error. 
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Lemma 4: The first two moments of Y(Z) can be expressed as 

o, n 
G(i) n(i),Ui. 



E[Y(0] = E [G(Z)R(0 + W(Z)] 



E [Y(Z)Y(Z) H ] = E [(G(Z)R(Z) + W(Z)) (G(Z)R(Z) + W(Z)) 



<r|,7r[b H (Z)] I + G(Z) b(Z)G H (Z) + 



where Tr [X] is the trace of matrix X and b(Z) — ^Sv + Sf + F(Z) F" (Z) 1. 
We make a Gaussian approximation on the multivariate observation vector to obtain: 
( \ 



Y(Z) ~ < 



CN Y(l);0,M44l + ^I 



/ 



CiV I Y(Z);G(i) IL(l),a 2 G Tr[b H (l)] I + G(Z) b(Z)G H (Z) - fi(l)fi H (l) + a^l | ,Hi. 

"mo" 



J Hl(!) 



The approximated distribution of Y matrix has the same structure as (20), and we can therefore utilise 
a similar procedure to obtain the decision rule: 

Lemma 5: Utilizing the results in Lemma 4 combined with Lemma 2 results in the LRT decision rule 
and Bayesian threshold as in Theorem 4. 

In making the GA, it is important to quantify the associated error with such a distributional assumption 
in evaluation of the LRT. Understanding the approximation error allows us to provide guidance on system 
design relating to the number of relay and the length of frames in order to mitigate errors in evaluating 
mis-detection and false alarms probabilities. 

Theorem 7: Under a Gaussian approximation to the distribution of the linearly transformed received 
signals Y(l), . . . , Y(M), where Y(m) = T (G^R^) , m = 1, • • • ,M, and T (•) is the linear 
standardization transformation, we obtain the Kolmogorov distance on all convex sets A £ A for 



c Y(l)+...+Y(M) . , 

M = Vm glven y 



400iV 1 / 4 E 



su PA& A\Pr (S M G A) - Pr (Z e A) \ < 



|Y(m) 



M 



where E 



|Y(m) 



r(f) 



, and Z ~ N (0, 1), 
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Proof: Using the result of Lemma 4, we transform the corresponding observation vectors Y(m) 
according to the following SVD decomposition of the covariance matrix of Y(m) given by S(m) = 
XJ(m)A(m)XJ(m) H . This produces the transformed i.i.d. random vectors given by 

Y(m) = U(m)A~ 1 / 2 (m) (Y(m) — E[Y(m)]). Having obtained i.i.d. vectors, we apply the multi-dimensional 
Berry-Essen bound [22]. To do this we only need to calculate E ||Y(m)|| 3 . Writing ||Y(m)|| = 

Y(!)(m)) 2 + . . . + (y^h) 2 , we have that (Y«(m)) 2 ~ xi (0) Vi E {1, . . . ,N}, and there- 
fore Y,n=i (Y( n )(m)) 2 ~ x% (0) and consequently, ||Y(m)|| ~ xn (0). The third moment of ||Y(m)||, 



which follows a xn distribution, is given by 2y2 '- . ■ 

\~2 ) 

Remark: considering convex sets of the form (— oo,x], Vx E R , the Berry-Esseen result shows the 
maximum error we can make under our Gaussian approximation of each observation vector and therefore 
provides a bound on the approximation error on the marginal likelihood used in the LRT 
Remark: the maximum error we can make under GA decreases at a rate of \[M (that is, the number of re- 



lays) for a fixed number of antennas, N. Furthermore, for a fixed number of relays, M, the approximation 
error becomes unbounded for increasing number of receive antennas, since ^ ^ oo. 

V 2 I 

B. Approximation of the Marginal Evidence via Laplace Approach 

In this section, a more accurate estimation of the marginal likelihoods than the GA is developed. This 
is based on the Laplace approximation [23]. The Laplace method can approximate integrals via a series 
expansion which uses local information about the integrand around its maximum. Therefore, it is most 
useful when the integrand is highly concentrated in this region. 

Under the full Bayesian paradigm, the evidence in (6) is obtained via the following marginalisation: 

p(y|?4) = y p(yl r >?4)p( r l?4) dr 

= ... j p{y\r ri ,H k ) ■ ■ ■p(y\r rM ,'Hk)p(ri\rik) ■ ■ ■ p {r M \H k ) dri . . .dr M , 

Jr\ Jr M 

where R = [R x (I),... , R M (l)] T - The densities in (34) can be expressed as: 

Y| (R = r; U k ) ~ CN ( Gr, (a 2 G ||r|| 2 + a^) l) (35) 



(34) 



CN(0,E V ), Ho 
R ~ F(r) ^ { (36) 

CN ( Fs,S F + Sy ) ,Hi. 
This integral is intractable and we shall approximate it via an application of the Laplace approximation. 
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To do so we begin by denning the following quantity (we discard the time dependency I here): 



h(r) = log(p(y|r)p(r)) . 



(37) 



This expression is now expanded using a Taylor series about its maximum a-posteriori (MAP) estimate, 
denoted by R = argmax r p(r|y). This is the point where the posterior density is maximised, i.e. the 
mode of the posterior distribution. Hence, we obtain 



h(r) =h R + r-R 



T dh R 



r-R 



T d 2 h R 



Q2 r 



r-R 



(38) 



dr 

(=0) at MAP location 

The second term in equation (38) cancels because at the maximum of h (r) (which is by definition what 
the MAP location represents), the first derivative is zero. 
Replacing h (r) by the truncated second-order Taylor series yields: 



h(r) &h(R\ +i (r-B.) H| 
where is the Hessian of the log posterior, evaluated at R: 



d 2 h(-R 



H 



d 2 lnp (r|y) 



drdr H 



r=R 



r=R 



We now concentrate on approximating the log of the integral in (34): 
logp(y)=log J p(y\r)p(r)dr 

= log J exp h( - r ^ dr 



Taylor seriei 



log|exp^ + K r - & ) TH ( r - a )dr 
h (R) + log f exp^^^dK 

<xCA r (R,H) 



I + ilog|2vrH| 



= logp (r) + logp (y|r) + |27rH 
Finally, the marginal likelihood estimate can be written as 

1 /2 

p (y) = v (?) p (y|r) |27rH _1 1 



-l|l/2 



(39) 



(40) 



(41) 



(42) 



The Laplace approximation to the marginal likelihood consists of a term for the data likelihood at the 
mode (second term of (42)), a penalty term from the prior (first term of (42)), and a volume term 
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calculated from the local curvature (third term of (42)). 

Under the Laplace approximation presented in (42), the LRT decision rule in (5) is approximated by 

MYvL) = nLmmii, f 7 (43) 

nti?(y(OI«o) «, 

where p{y(l)\T-Lk) is the Laplace marginal likelihood approximation under the fc-th hypothesis. The major 
difficulty in evaluating (42) is the requirement to evaluate the MAP estimate R under each hypotheses. 
This task is nontrivial as it involves a non-convex and non-linear optimisation problem. We derive the 
MAP estimate for this scenario via the Bayesian Expectation Maximasation (BEM) methodology, see the 
derivation in the Appendix . 

VI. Simulation Results 
In this section, we present the performance of the proposed algorithms via Monte Carlo simulations. 



A. Simulation Set-up 

The simulation settings for all the simulations are as follows: 

• The prior distribution for all the channels is Rayleigh fading, and the channels are assumed to be 
both spatially and temporally independent. 

• We define the receive SNR as the ratio of the average received signal power to the average noise 
power 

1 



SNR = 10 log 



Tr 


K 


"(G(I)F(IMI)) (G(Z)F(0 S (0)^ 






Tr 


E 


"(G(l)V(O + W(Z)) (G(Z)V(O + W(l)) H ] 



10 log 



°V + M °W 



The SNR is set to dB. 

The results are obtained from simulations over 100, 000 channels and noise realisations for a given 
set of N, M and L. 

For the Laguerre series expansion, the order of the series expansion was set to p = 100. 



B. Study of detection probability Vs. frame length 

In this section we study the relationship between the ability to detect the presence of a signal in a 
spectrum sensing problem as a function of the length of the frame, L. We undertake this study in two 
different scenarios, the first involves perfect CSI according to Section III and the second involves partial 
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CSI according to Section IV. We set the channels uncertainty F(7) = and <7 F = 1, thus only prior 
information is available for the F channels. 

In presenting results we fix the false alarm rate pj to 10%. We repeat this study for a range of values of 
the number of receive antennas, N € {1, 2, 4, 8}. The results are depicted in Fig. 2 and they demonstrate 
the following key points: 

1) for all frame lengths, as the number of receive antennas is increased, the probability of detection 
improves as expected; 

2) for all frame length, the detection probability under perfect CSI always outperforms significantly 
the performance of the model with partial CSI; 

3) asymptotically in the frame length, L, the probability of detection for any number of receive 
antennas converges to 1, with different rates, depending on ./V; 

4) such a study provides generic performance specifications that allow us to obtain the same detection 
probability for different combinations of frame length and number of receive antennas. For example, 
with L = 10 and N = 1, this will be equivalent to L = 3 and N = 4. 

5) it also guides system design that for a given desired probability of detection, we see the saturation 
point, after which, increasing the frame length delivers negligible improvement. 

Evaluation of the Gaussian approximation 

The accuracy of the GA in Section V-A is bounded via a multi-dimensional Berry-Esseen inequality 
in Theorem 7. Here we study this accuracy using a graphical Q-Q plots of each element of Y(7) as a 
function of the number of relays M. The results are presented in Fig. 3 and demonstrate that for a fixed 
frame length and number of receive antennas, as one increases the number of relays M, the Gaussian 
approximation that we made in Lemma 4 improves. We see that in the setting of partial CSI which is 
relevant to practical scenarios, the number of relays required before one can make a reasonable Gaussian 
approximation is around 8. 

C. Comparison of Detection Probability under different LRT Statistic Approximations 

In this section we present a comprehensive comparison of the distributional estimators derived for 
the LRT test statistic in order to evaluate the probability of detection. This is undertaken in a range of 
different scenarios and we compare the distributional estimates under different levels of CSI versus the 
best case scenario bounds. The comparison is undertaken between: 

1) the analytic evaluations of the probabilities of detection and false alarm under the setting of perfect 
CSI, according to results obtained in Theorem 2 (denoted by: CSI Theory); 
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2) the Monte Carlo based empirical estimation of the probabilities of detection and false alarm under 
the setting of perfect PBS-relays CSI, and perfect relays-SBS CSI, according to the decision rule 
derived in Theorem 1 (denoted by: CSI empirical); 

3) the analytic evaluations of the probabilities of detection and false alarm under the setting of 
imperfect PBS-relays CSI, and perfect relays-SBS CSI, according to the Laguerre series expansion 
density approximations derived in Theorem 5 in (22a-22b) (denoted by: P-CSI Laguerre); 

4) the Monte Carlo based empirical estimation of the probabilities of detection and false alarm under 
the setting of imperfect PBS-relays CSI, and perfect relays-SBS CSI, according to the decision rule 
derived in Theorem 4 (denoted by: P-CSI empirical); 

5) the Monte Carlo based Gaussian approximation of the probabilities of detection and false alarm 
under the setting of imperfect PBS-relays CSI, and imperfect relays-SBS CSI, according to Lemma 
4 applied to the decision rule derived in Theorem 4 (denoted by: PP-CSI Gaussian); 

6) the Monte Carlo based Laplace approximation of the probabilities of detection and false alarm under 
the setting of imperfect PBS-relays CSI, and imperfect relays-SBS CSI, corresponding decision rule 
also derived (denoted by: PP-CSI Laplace). 

The scenarios we consider involve varying the number of receive antennas N and the number of relays 
M, for a fixed frame length L = 1 and a fixed SNR of dB. The Receiver Operating Characteristic 
(ROC) curves are presented in Figs. 4- 9, for each of these comparisons. The following summary details 
the key points of this analysis: 

1) In all study combinations of N and M, the probability of detection for each probability of false 
alarm, had an ordering of algorithmic performance, in agreement with theory, given by: 

i. Optimal performance under perfect CSI. This results in the theoretical upper bound of Theorem 
2 which agreed exactly with the Monte Carlo estimate under this scenario. 

ii. This was followed by the results of the imperfect PBS-relays CSI, and perfect relays-SBS 
CSI which were obtained under the Laguerre approximation and again compared to a Monte 
Carlo simulation estimated. 

iii. Finally the results of the approximations when least information is known, imperfect PBS- 
relays CSI, and imperfect relays-SBS CSI which were obtained under the Laplace approxima- 
tion and the Gaussian approximation. The Laplace approximation outperformed the Gaussian 
approximation in situations in which the distribution of the test statistic was not close to 
Gaussian. 
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2) In all the examples the Laplace approximation outperformed the Gaussian approximation or was 
directly comparable in performance as the Central Limit Theorem became viable, i.e. when M was 
large, as presented in Fig. 6. 

VII. Conclusions and Future Work 

In this paper we developed a framework for spectrum sensing in cooperative amplify-and-forward 
cognitive radio networks. We developed the Bayesian optimal decision rule under various scenarios of 
CSI varying from perfect to imperfect CSI. We designed two algorithms to approximate the marginal 
likelihood, and obtained the decision rule. We utilised a Laguerre series expansion to approximate the 
distribution of the test statistic in cases where its distribution can not be derived exactly. Future research 
will include comparison of the Laplace method to other low complexity approaches, such as the Akaike 
and Bayesian information criteria. 
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Appendix 

Proof of Theorem 6 

Proof: Consider the normalised product of two independent normally distributed random variables 
defined as 



(JxCTY 

where X ~ N (X, a\) and Y ~ N (F, cr^), and define p x = £ and p Y = 
The expectation of Z is given by 

- . n E\X}E\Y] 
Z = E[Z]= l _ J _ 1 J = pxPY, 



(44) 



(JXCTY 



The variance of Z is given by 



a\ = E [Z 2 ] - E [Zf 



E [X 2 ] E [Y 2 



z z 



a 2 x +X z ) (a 2 v + Y 2 



•y 



°\°Y 



E[Z] 2 =l + p 2 X +Py. 



We define Z = (Z -~Z) /a z and derive the MGF of Z using (31): 



M~ z (t) =M {z _- 2)laz (t) = E 



exp 



<?Z 



-—t ; 



, z 

exp < — t 



exp 



exp 



J \CTZ 



PxpY 



yJl+Px+pl 



I 2 , 2 \,2 



2pxPY f 



t > exp 



2 1 



I+Px+Py 



2 PXPY t 



exp 



pxpvt 



(45) 



2 1 



V 1+ p 2 x + pI 



l+p x +p Y 



{p 2 x+p Y )t 



2 PX PY 1 



2 PXPY t 



exp 



2 1 



f 2 , 2 \.2 



2 PX PY 1 



exp 



2 1 
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Finally, we define a = , 1 and take the limits px — > oo or py — > oo to obtain the following 

V x +Px+Py 

standard normal distribution: 

eX P 1 2(l-a 2 ) I ( t 2 } 

lim Ms = lim ^ 1 ; } - = exp { — }, (46) 

Px -> oo 

/9y — )• OO 

which is the MGF of iV (0, 1). Therefore, we obtain that Z ~ (Z, cr|). ■ 
Deriving the MAP estimate of R in (42) 

The MAP optimisation problem can be written as 

R = arg maxp (r|y) = argmaxp (y|r) p (r) , (47) 

r r 

where p (y|r) and p (r) are defined in (35)-(36). Then, the MAP estimate is the solution for the following 
optimisation problem, where for simplicity we remove the time dependence I: 

R = arg max p (y|r) p (r) 



r 



lk-°HI 2 > 1 ( Ik-* 



argmax^ ^expV "cii'i^+vy x — — g-expv " R / , (48) 



where 



p(yk) 



n : R = 0, a^ = al 

U x : R = Fs,<xf l = 4 + ^v- 



p(r) 



(49) 



Problem (48) is non-linear and non-convex. We shall utilise the Bayesian Expectation Maximisation 
(BEM) methodology to solve it efficiently under each hypothesis. The BEM algorithm (see [24], [25]) 
is an iterative method that alternates between an E step, which infers posterior distributions over hidden 
variables given a current parameter setting, and an M step, which maximises p (y, G, R) with respect to 
R given the statistics gathered from the E step. The BEM can be easily evaluated using the following 
iterative steps, at iteration (n + 1): 

E Step: L (R) = E G|y;ft „ [logp (y, G, R)] (50a) 
M Step: R" +1 = arg max L (R) (50b) 

R 
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The E Step can be expressed as: 
L(R)=E G|Y . a „ [logp(Y,G,R)] 

= E G|Y;R" P°gP ( Y I G ' R ) + l °SP (R-PK)] 
= E 



J G|Y;R" 



(7 



-J-||Y-GR|| 2 
w 



(7 



— ||r — r|| 

R 



+ constant 



= — (2Y H RE G|Y . fi „ [G] - R T E G|Y . fi „ [G T G] r) - \ (R^R - 2R H R) + constant. 

(51) 

where constant contains all terms that are independent of R. 

The conditional expectations in (51)can be evaluated using Bayesian MMSE as follows: we first re-write 
the observation model (3) as: 

Y = GR + W = (R T ® I) vec [G] + W = QT + W, (52) 



where we define Q 



R 



T 



I), r = vec [G], and is the Kronecker product, and vec [•] is the vector 



obtained by stacking the columns of a matrix one over the other. Since Y and T are jointly Gaussian, 
the Linear MMSE is also the MMSE estimator (see [26]). The LMMSE can be expressed as 



E r|Y-R- = e [r] + e [rY ;/ ] e - 1 [yy"] (y - e [y]) 



r + E [TY H ] E" 1 [YY H ] (Y - SIT) 



- n H (y - nT) 

T + 



+ -T 1 



where r = E [T]. Next, we evaluate the covariance matrix: 

Cov r|Y-R» t r l =E[rr H ]-E [r y h ] e- 1 [y y h ] e [y t h ] 



By rearranging the above expressions, we obtain 
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R n ) =E G|Y;ft „ [G] = G 
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R n 





Y - GR" R" 
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<J? 2 (Y, R") 4 E G|Y . ft „ [G H G] = $! (Y, R") (Y, R n ) + a 2 G N 
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(55b) 
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Using (55a-55b), (51) can be expressed as 

L(R) = \- (2Y^R$i (Y, R n ) - R T $ 2 (y, R") r) - — (R H R - 2R H R) + constant. (56) 

The M Step is obtained by setting to the derivative of L(R) with respect to R: 

R n+1 = arg max L(R) = ($ 2 (y, R fc ) + ^-ij (V (y, R fc ) T Y + • ( 57 ) 

The BEM algorithm requires that R™ +1 is initialised at n = 0. The simplest option is to initialise it to 
the prior, that is R° = R. 



Case 


PBS-relays (F) 


relays - SBS (G) 


Section 


Decision rule 


Performance analysis 


I 


V 


V 


III 


Exact 


Exact analytic 


II 


X 


V 


IV 


Exact 


Analytic approximation 












via Generalized Laguerre polynomials 


III 


V 


X 


V 


Special case of IV 


see Section IV 


IV 


X 


X 


V 


Analytic approximation 


Simulation 










via Laplace integrals 




IV 


X 


X 


V 


Analytic approximation 


Simulation 










via Moments matching 





TABLE I 

Summary of proposed solutions based on CSI knowledge 



Primary User 
Base Station 



Fm*. 




Secondary User 
Base Station 



Relay M 



Fig. 1. System model of Cooperative Cognitive Radio network with M relays and a multiple antenna receiver 
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Fig. 2. Probability of detection at pf — 0.1 for the cases of perfect CSI (Section III) and imperfect CSI (Section IV) 




4-3-2-101234 
Slandard Normal Ouantiles 



M-i 




Slandard Normal Ouantiles 



11=100 




4-3-2-101234 
Standard Normal Qoantiles 



Fig. 3. Q-Q plot of the normal approximation per Section V-A for different number of relays 
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Fig. 7. Probability of detection vs. probability of false alarm for N = 8, M = 2, L = 1 
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Fig. 9. Probability of detection vs. probability of false alarm for N = 8, M = 8, L = 1 
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